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SUMMARY 


In recent years the surface-source method of calculating potential flow 
about arbitrary bodies has been developed extensively and has proved to be a 
useful tool in a wide variety of low-speed design applications ranging from 
simple shapes to complicated Inlets with centerbodies, multielement airfoils, 
and wing-fuselage-pylon-nacelle combinations. Two-dimensional, ax 1 symmetric, 
and three-dimensional methods have been developed. While the method is gener- 
ally quite satisfactory, increases in computational speed and accuracy are 
desirable for certain applications, particularly interior flows and exterior 
flows about complicated multiple-body combinations. Such improvements can be 
realized by refining the formulation. In the basic method the profile curve 
of a two-dimensional or axisymmetric body is approximated by a large number of 
straight-line elements over each of which the source density is constant. The 
so-called higher-order refinement consists of using curved surface elements 
and a source density that varies over an element. This report describes the 
analysis for the axisymmetric case where the next-simplest approximation is 
used. Specifically, the surface elements are parabolas, and the source density 
varies linearly in arc length over each element. Calculated results for a 
series of test cases are presented and compared. For both exterior and interior 
flows the higher-order formulation yields very significant decreases in comput- 
ing time and increased computational accuracy. 
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LIST OF SYMBOLS 


radius of a ring source or vortex 
axial location of a ring source or vortex 
half the curvature of an element 
cy 

integer subscripts denoting particular elements 

total number of elements on a body 

arc length along an element 

half the total arc length of an element 

s'y 

surface velocity 
free-stream velocity 
velocity components 

velocity due to a ring source or vortex 
perturbation velocity at i-th control point 

velocity at i-th control point due to a unit value of j-th source 
density with all other values zero 

velocity at i-th control point due to a parabolic source density 
distribution on the j-th element 

axial coordinate 

coordinates of a control point 

radial coordinate 

slope angle of the tangent to an element at the control point with 
respect to the x-axis 

length of the chord of an element 

perpendicular distance from control point to chord of' an element 



8 

y 

n 


a (1 > 

, (2) 

♦ 


circumferential angle about the x-axis 
vorticity strength 

Cartesian coordinates with origin at the control point of an element 
and directions, respectively, parallel and perpendicular to the 
tangent line 

source density 

derivative of source density with respect to arc length 

half the second derivative of source density with respect to arc 
length 

velocity potential 
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INTRODUCTION 


Recent years have seen the development of very general surface singularity 
methods for the calculation of potential flow about arbitrary configurations [1], 
Moreover, these methods have been applied successfully to a large number of 
practical design problems of low-speed flow [1], [2], The most common such 
method is the so-called surface-source method [2], which utilizes a source 
density distribution over the surface of the body about which flow is to be 
computed. Application of the boundary condition of zero (or prescribed) normal 
velocity on the body surface theoretically yields a Fredholm integral equation 
of the second kind for the source density. Once the source density is known, 
all other quantities, such as flow velocities and pressures, may be obtained by 
integration. Separate procedures have been developed for calculating flow about 
two-dimensional bodies, axisymmetric bodies, and three-dimensional bodies, 
respectively. For the case of axisymmetric bodies the flow itself does not 
necessarily have ‘ ~''mmetric, but it may be a case of cross flow for 

which the free-stream direction is perpendicular to the body's symmetry axis 
or a case of rotation of the body about an axis perpendicular to its symmetry 
axis [2]. 

To impleffient this method for the computer, various approximations must be 
made. In particular, both the body shape and the source density distribution 
must be approximated in a form suitable for machine computation. During most 
of the development of the surface-source method the profile curve defining an 
axisymmetric (or two-dimensional) body has been approximated by a large number 
N of small straight-line elements, which form an inscribed polygon. Moreover, 
the source density has been assumed to be constant over each straight-line 
element, although it varies from one element to another. This reduces the 
problem of determining the source density distribution to that of determining 
a finite number of values of the source density — one for each element. One 
point of each element, the midpoint for a straight-line element, is selected 
as the control point where the normal -velocity boundary condition is to be 
applied. Formulas have been derived that give the velocity at any point due 
to a unit value of source density on a straight-line element. From these 
formulas a matrix of vector velocities induced by the elements at the control 
points can be obtained. Then the integral equation is replaced by a set of 
linear algebraic equations forthe values of source density on the elements. 

The coefficient matrix of this set of equations consists of the set of normal 
velocities induced by the elements at each others' control points, which is 
obtained by taking normal components of the basic induced velocity matrix. 
Finally, surface velocities at the control points are obtained by a matrix 
multiplication of the tangential components of the induced velocity matrix 
with the column of values of source density. The two main parts of the compu- 
tation are the calculation of the velocities that comprise the induced 
velocity matrix and the solution of the linear equations for the values of 
source density. If a direct elimination solution is used, the computational 
magnitude of solving the linear equations is proportional to N^. 

The above described procedure, which uses flat surface elements and a 
piecewise-constant source density distribution, is designated the base method. 

It has proved satisfactory in a wide variety of design applications [1],[2]« 
However, it is evident that more elaborate procedures can be formulated and 
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that these should give a higher accuracy for a given element number N and 
thus an equal accuracy with a smaller element number. Moreover, if it is 
properly implemented, a more elaborate procedure would require only a slightly 
greater computing time than the base method for a given N. Because of the 
rapid variation of computing t‘ime with N, a reduction of computing time for 
a given accuracy should be possible. This has been successfully accomplished 
for the case of two-dimensional bodies [1], [3], [4], [5], in which case the 
application of main interest is a multielement airfoil. Here the case of 
axisymmetric bodies is considered, where the application of main Interest is 
an inlet, possibly with centerbody and ring vanes. The technique employed is 
basically the one used in two dimensions [3], which employs curved surface 
elements and a source density that varies over the element. Such an approach 
is designated a higher-order implementation. For axisymmetric bodies these 
variations refer to the profile curve defining the body. The circumferential 
variations around the symmetry axis of element geometry and source density 
are similar to those of the base method [2]. Namely, a complete surface ele- 
ment is a portion of a conoid (a cone frustum in the flat-element case), and 
the source density is independent of circumferential angle for axisymmetric 
flow and is proportional to the cosine of the circumferential angle for cross 
flow and for rotation. 

Use of higher-order implementations brings up the question of the consist- 
ency of the orders of the approximation used for the element geometry and the 
approximation used for the source density. It is known [1], [3], that consist- 
ency is obtained when the polynomial expressing the element shape has a degree 
one higher than that defining the source density. Thus, consistent formulations 
include: (1) flat-element constant-source, (2) parabolic-element linear-source, 

and (3) cubic-element parabolic-source. In this paper parabolic elements are 
used together with a piecewise-parabolic source density. The extra inconsistent 
parabolic term in the source density has been included to determine whether or 
not it gives increased accuracy in low-curvature regions. It turns out that, 
as the theory predicts, little or no gain in accuracy is obtained by including 
the parabolic source term. 
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SURFACE ELEMENT GEOMETRY 


The profile curve of the body about which flow is to be computed is 
specified as a table of coordinates for N + 1 points (xp y-j), each of 
which is presumably exactly on the contour. By this means the contour is 
divided into N elementary arcs as shown in Figure K On each arc a control 
point is selected by the following criterion: the normal projections of the 

endpoints of the element, (xp y-j) and (xpp yp^), on the line tangent 
to the arc at the control point are equidistant from the control point. The 
slope of the tangent line at the control point is defined as the slope of the 
element. If this tangent line is taken as the horizontal axis of a in- 
coordinate system with the control point as origin, the elementary arc is as 
shown in Figure 2, The equation of this arc may be written as a power series 

2 

n = c£ + . . . j-j 

The arc length s along the arc measured from the control point is given by 


ds = [Vl + (2 c 5) 2 + ...]ds (2) 

= [1 + 2c Z E 2 + .]d£ 

Thus 

2 2 3 

s = % + f cV + ... (3) 

The basic reference coordinate system in which the body is defined has its 
x-axis as the symmetry axis of the body„ Let a point of the element have 
coordinates x = b y = a in this system then 

— 2 

b = X + (cosa)c - c(sina)s + ... 

- 2 ' ' 

a = y + (sina)£ + c(cosa)c + ••• 

where a is_the^ slope angle of the £~axis (tangent line) with respect to the 
x-axis and x, y are reference coordinates of the control point. 

In principle, approximations of arbitrarily high order could be generated 
by retaining surfficient terms in the above series. For present purposes the 
element is assumed to be parabolic, and the above series are terminated with 
the terms shown. For this approximation the control point lies on the perpen- 
dicular bisector of the straight line between (xp y-f) and (xpp Yi+l) 
and the slope of the element equals the slope of this straight line. Also, by 
inspection of (3) it can be seen that to this order of approximation s may be 
substituted for £ in (4). A circle is passed through the points (xpp } > 

(xp yih and (x.^, y. + -j), and another circle is passed through the points 
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(xj, yj), tx. +1 , y 1+1 ) and (x i+ 2 » yi+?)* The half-curvature c of the para- 
bolic element, which enters equations (T) and (4), is set equal to half the 
geometric mean of the curvatures of these two circles. Requiring the parabola 
to pass through the endpoints (x-j , yj) and (x^i , y-j+i ) then uniquely 
defines the parabola and establishes the coordinates x, y of the control point. 
For the first and last elements of a body one circle is not defined and c is 
set equal to half the curvature of the remaining circle. 

Explicit formulas for the above procedure are contained in Appendix A. 
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INDUCED VELOCITY MATRICES 

The basic calc lational task of the fl (^calculation method is to compute 
the How velocities induced by the elements at each others' control points. 

Ther » are three tyf :s of induced velocities corresponding to three types of 
surface singularit (1) a constant-strength ring source, which is appropriate 
for . ixi symmetric 1 w; (2) a ring source whose strength is proportional to the 
cosi le of the circ inferential angle, which Is appropriate for cross flow and 
for jody rotation; and (3) a constant-strength ring vortex, which is used to 
produce certain auxiliary solutions appropriate to axl symmetric flow about 
ring -airfoils and inlets. The velocities in the induced velocity matrices are 
obtained by integrating the ring source or vortex formulas over an element. 
Specifically? let a ring singularity (source or vortex) have a radius a and 
lie in the plane x = b with its center on the x-axis. Then the velocity at 

the point (x,y) due to this singularity is 

7 = y, b, a] (5) 

Specific formulas for the two types of ring source are contained in [2] and 

derived in [6], [7]. Many of the functions appearing in the expressions are 

repeated, particularly the specific complete elliptic integrals. The velocity 
components of a constant-strength ring vortex are related to those of a constant 
strength ring source by [8] 


V x (vortex) = V x (source) + ^ V^(source) 

( 6 

Vy(vortex) ■ V x ( source) - V^(source) 

Let the source density on the j-th element be denoted a.(s). The velocity 
induced by this source density at control point of the i*th element (x-, y.) 
is obtained by Integrating over the element, that is 


V*. = /vtV y v bj(s), a j (s)]aj(s)ds (7) 

AS. 

J 

where as. is the total arc length of the j-th element. Equation (7) applies 
to both types of source singularity. This equation also applies to the vor- 
ticity singularity if is replaced by the vorticity strength pj(s). 

In evaluating the integral of (7), b and a are given by (4) with s replac- 
ing 


The source density may be written as a power series 


a.(s) = a. + oP's + a^S 2 + ... 
J J J 


(8) 
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where crj, <jP^ and oP are Independent of s and represent, respectively, 
the value, first derivative, and half the second derivative of the source 
density at the control point of the j~th element. In the present analysis the 
source density is assumed to be at most parabolic, and (8) is terminated with 
the terms shown. Thus (7) becomes 

% = °j /ft*!’ V b j(s), a.(s)]ds 

ASj 

+ cjP) J Vt7., y^, bj(s), aj(s)]sds (9) 


or 




b j ( s ) , aj(s)]s 2 ds 


+ V< 2 U 2 > 

10 1J 0 J 1J J 


(10) 


The vorticity distribution y - ( s ) is taken as constant so that only the first 
terms of (9) and (10) are present. In fact, the vorticity is handled exactly 
as in the base method [1], except that the integral is over a curved element 
rather than a flat one. The form of equation (10) makes it easy to investigate 
the effectiveness of retaining various terms in the source density expansion. 
The relative effectiveness of flat and curved elements can be investigated by 
setting c either zero or nonzero in equation (4). 

For j f i the integrals in (9) are evaluated by numerical integration 
using Simpson's rule with a variable number of ordinates [2], It is evident 
from (9) that the three integrands are very similar and can be conveniently 
calculated together. For j - i the basic ring-source (or vortex) velocity 
Vtx,y,b,a] has a singularity of order 1/s at s = 0. Accordingly, resort 
must be made to a series expansion in powers of s so that the singularity 
may be cancelled analytically. These expansions are similar to those used in 
the base method [2]. The quantities vtx,y,b,a]s and V[x,y,b,a]s2 are not 
singular but for consistency they are evaluated by series expansions for 
j = i. Derivation of these expressions is straightforward but tedious and is 
not pursued here. Formulas for these expansions are contained in Appendix B. 
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ORGANIZATION OF THE CALCULATION 


In the most general case all three integrals of equation (9) must be 
evaluated for both the axi symmetric and the cross- flow type of source density. 
Thus including the vorticity, which requires only the first integral of (9), 
there are seven vector integrals in the higher-order analysis as opposed to 
three in the base method. In axi symmetric flow there are two velocity compon- 
ents, axial and radial, while cross flow has an additional circumferential 
component. Thus the total number of scalar integrals can be seventeen rather 
than the seven of the base method. However, the final number of scalar N x N 
induced velocity component matrices that must be stored and used is the same in 
the higher order and base methods, namely five. 


As described above, the normal velocity boundary condition is applied at 
the control point of each element to produce a number of linear equations equal 
to the number of elements. However, the variation of source density over an 
element is described by three parameters, oj, oj^), and 
tives aj(l) and aj(2) are expressed in terms of values of J 
entiating the parabola through the three values o. , a., 

, . J ^ I J 

( 1 ) 


and 


The deriva- 
by differ- 
ed . That is 
J+l 


Vm + E j°j + F j°j+i 


OD 


( 2 ) . 


G -a • i 
J J-l 


Vo 


Vj+I 


where the coefficients in (11) are standard numerical differentiation formulas 
and depend only on the lengths of the three elements [3]. For the first (or 
last) element of a body the parabola that is differentiated passes through the 
first (or last) three values of a, and the formulas of (11) are modified 
accordingly. This last feature could introduce error if a smooth contour is 
defined with the first and last points coincident, as for example, a torous. 
However, this case is too rare to be of great significance. The above pro- 
cedure Is appropriate for all inlets and ducts and also for simply-connected 
bodies, which are input from the upstream stagnation point to the downstream 
stagnation point. 


From (10) and (11) it is evident that the velocity induced at a point by 
an element is a linear combination of three neighboring values of a. The 
matrix that is needed for subsequent calculations is the one giving the veloc- 
ity at each control point due to each value of source density. As the calcu- 
lation proceeds to calculate the velocity at a control point due to successive 
elements, the velocity induced by the j-th element is not associated entirely 
with o- aj in the base method, but certain portions are associated with 
oj-1 » cn , and g. + , in the obvious way. When all elements have been accounted 
for, the velocity J at a control point due to each value of source density has 
been formed as the sum of contributions from three elements except for those 
due to the first and last values, to which two elements contribute and then 
those due to the third and third-to-last values, to which four elements contri- 
bute. The result of this phase of the computation is a matrix V^j, such that 
the velocities at the control points due to the body are 
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1 = 12 , 


• •• » 


N (12) 


N 

v. = > V. .a • , 

1 L~i 13 j 

There is a single vector matrix Vy for the axl symmetric source densities, 
one for the cross-flow source densities, and one for the axi symmetric vortic- 
itles, just as In the base method* Moreover, as in the base method, the vortex 
velocities are not saved individually. Instead the velocities produced by all 
elements at a control point are added together to give the velocity at that 
control point due to a vortex sheet of constant unit strength. This auxiliary 
onset flow is used in certain inlet and ring-wing applications. 

Thus the result of this phase of the calculation is a set of matrices • 

and vorticity onset-flows equivalent to those of the base method. In fact, J 

all subsequent calculations [2] are identical for both base and higher-order 
methods and do not depend^on how the V-jj matrix was produced. In particular, 
the normal component of V-jj Is the coefficient matrix of the linear equations 
for the values of source density* The right sides of these equations are the 
negatives of the normal components of the onset flows, either uniform or other- 
wise. Once solutions have been obtained, the flow velocity at each control 
point (or any other point) is calculated for each onset flow by adding a sum of 
the form (12) to the onset- flow velocity at that point. 

Explicit coefficients for the numerical differentiation formula (11) are 
presented in Appendix C. 
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DISCUSSION OF RESULTS 

Comparison of Calculated Results with Analytic Solutions 

To determine the effectiveness of the higher-order technique and to evalu- 
ate the Importance of the various terms in the expansion, a considerable number 
of calculations have been performed for a sphere and for a 8-to-l prolate spher- 
oid, for which analytic solutions are available. The first conclusion that can 
be drawn from these calculations is that the addition of the quadratic source 
density term a A 2 ) never yields an appreciable increase In accuracy regardless 
of whether flat J or curved surface elements are used. In every case the solution 
that utilizes only aj, and cr-j(l) is virtually Indistinguishable from the one 
that also includes In what follows, the solution that utilizes curved 

surface elements and a linearly varying source density is denoted "higher order" 
while the flat-element constant-source solution is denoted the base method. 
Solutions obtained with other combinations of- terms are labeled explicitly. 

Calculations were performed for a sphere represented by 60 equal-length 
elements of 3° subtended angle. Four solutions were obtained. In addition to 
the higher order and the base method, two inconsistent solutions were obtained: 
that using the curved-elements constant-sources and that using flat-elements 
linear-sources. The higher-order solution was also calculated for a 12-element 
sphere whose elements each subtend a 15° angle. Results for a uniform onset 
flow parallel to the x-axis are presented in Figure 3, which shows differences 
between calculated and analytic surface velocities. The importance of mathe- 
matical consistency is evident. Accounting for either source variation or 
element curvature separately produces no improvement on the base method. How- 
ever, when both are accounted for, as in the higher-order solution, the result 
is a large gain in accuracy — about two orders of magnitude. Even the 12-element 
higher-order solution is an order of magnitude more accurate than the 60-element 
base method and also requires an order of magnitude less computing time. 

If the uniform onset flow is parallel to the y-axis, i.e., a "cross flow," 
the calculated results are slightly different from those for the case of onset 
flow parallel to the x-axis, because in the former case the conoidal surface 
elements do not have the same axis of symmetry as the flow field. Errors in 
calculated surface velocities are shown in Figure 4. Figure 4a shows velocities 
along the profile curve in the Diane containing the onset flow vector and the 
body's symmetry axis (xy-plane). Figure 4b shows velocities along the curve 
in the plane perpendicular to the onset flow and containing the body's symmetry 
axis (xz-plane). In this last plane the analytic velocity has a constant value 
of 1.5 and is parallel to the onset flow vector. The gain in accuracy from use 
of the higher-order solution is not as impressive in Figure 4a as it is in Fig- 
ure 3, but it is still present. In particular the 12-element higher-order 
solution is still at least as accurate as the 60-element base method. The 
ineffectiveness of the two inconsistent solutions is quite clear. The results 
of Figure 4b are similar, with one exception, to those of Figure 3. Both the 
60-element and the 12-element higher-order solutions are much more accurate 
than the 60-element base method. However, unlike the case of Figure 3, the 
gain in accuracy of Figure 4b is due largely to the use of curved elements 
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while the use of a variable source density has virtually no effect. However, 
the fact that this is true for one velocity component out of three does not 
change the basic conclusion that In general a consistent formulation is required 
to obtain a significant increase in accuracy. 

The relative effectiveness of the two methods of solution for an 8-to-l 
prolate spheroid Is illustrated In comparing results for the base method using 
60 elements with results for the higher-order method using 30 elements. For 
these element numbers the higher-order method is approximately four times 
faster than the base method. Somewhat different distributions of elements are 
used for the two methods. The distribution used for the higher-order method 
concentrates elements in the high-curvature region to a greater degree than that 
for the base method. Correspondingly, the higher-order solution uses very 
few elements in the low-curvature region, but it is able to maintain accuracy 
in this region because of the use of variable source density. Generally, each 
method is used with the element distribution for which the best solution is- 
obtained, (Similar experience is reported in [3].) Differences between cal- 
culated and analytic surface velocities are shown* in Figure 5 for the axi sym- 
metric case where the onset flow is parallel to the body's symmetry axis (x-axis) 
and in Figure 6 for the cross-flow case where the onset flow is perpendicular 
to the body's syronetry axis (parallel to the y-axis). To put these results in 
perspective, the maximum value of surface velocity is about 1.0293 times free- 
stream velocity for the axisymmetric case of Figure 5 and about 1.9447 times 
free-stream velocity for the cross-flow case of Figure 6, In the xz-plane. 
Figure 6b, the velocity is parallel to free-stream and has a constant magnitude 
of 1.9447. The curves of Figures 5 and 6 egiphatlcal ly show that in addition to 
being much faster the higher-order technique is also much more accurate than 
the base method. The substantial gain in accuracy that can be achieved by use 
of the higher-order method for cases of smooth convex bodies is in marked con- 
trast to the two-dimensional case [3], 

Interior Flow in Ducts 


In two dimensions [4] the greatest gains in accuracy from use of the higher- 
order method occur for interior flow in ducts. To investigate the behavior of 
the present axisyrrmetric case, the first geometry selected is the analog of one 
previously considered in [3], namely the case of a uniform onset flow into the 
closed duct shown in Figure 7 (an interior hemisphere cylinder),, The flow 
inside the duct should be virtually stagnant, and the average axial velocity 
component at any axial location should be exactly zero. Thus, the average 
axial velocity is a measure of the calculational error 0 The base method with 
55 elements gives a typical error of 15 percent of free-stream velocity. Use 
of the higher-order technique reduces this error by a full two orders of 
magnitude. 

A more practical case is that of a contracting duct of area ratio 16 as 
shown in Figure 8. The contracting section is a portion of a sine wave. By 
continuity the net flux of fluid at every axial station should be identical, 
so changes in this flux represent calculational error. Net fluxes have been 
calculated at five stations as shown in Figure 8 and normalized with respect 
to the flux in the large constant-diameter section D It is evident that the 
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use of the higher-order solution reduces the maximum error in flux by a factor 
of about 25. 


Inlets 


The most frequently-occurring type of axl symmetric body is an inlet. 
Simulating an arbitrarily prescribed mass flow ratio through the inlet requires 
use of an artifice that is described in detail In Appendix D. Basically, three 
flow solutions are required: that about the empty Inlet in a uniform onset flow 

at zero angle of attack (Figure 9a ) s the pure cross-flow due to a uniform stream 
at 90° angle of attack (Figure 9b), and a static solution, where there is flow 
through the inlet but no onset flow. This last solution can be obtained either 
by the use of an interior suction surface (Figure 10a) or by use of surface 
vorticity (Figure 10b) . The first two solutions are no more difficult than those 
for flow about simple closed bodies. It is the static solution that can lead to 
numerical difficulty and may require large element numbers. 

Several inlet goemetries were studied with both the higher-order solution 
and the base method for two or three different element numbers , Since results 
for all inlets are similar, discussion here is concentrated on the results for 
the inlet-centerbody combination shown in Figure 11. The control station where 
mass flow ratios are evaluated is at x = 13.4 and the inlet is terminated at 
x = 44, about two diameters. Figure 12a shows an input point distribution 
typical of those customarily used to define such an inlet - a total of 244 
elements. The element number has been reduced by deleting every other input 
point in the forward region, but the same element size is maintained on the 
afterbody. The resulting input point distribution, which corresponds to 149 
elements, is shown in Figure 12b. The same procedure was applied to the dis- 
tribution of Figure 12b to produce the 103-element distribution of Figure 12c s 
which has only one-fourth the density in the nose region as that of Figure 12a. 

The most surprising result of the calculations for all inlets and for both 
the higher-order and the base method is that the calculated velocity distribu- 
tions do not change much with element number. Thus, if a user contemplated a 
series of cases, it would be profitable for him to do a little initial experi- 
mentation with element number and select the lowest possible. 

With regard to flow continuity inside the inlet, the higher-order solution 
is an order of magnitude more accurate than the base method, but the latter is 
accurate enough for most purposes. The only real advantage in this regard for 
the higher-order solution is that its control station may be located anywhere, 
while the base method should have its control station as far forward as possible. 

The higher-order solution and the base method give virtually identical 
results for the two solutions of Figure 9. This was not unexpected. However, 
for the static solution of Figure 10, the higher-order solution offers an 
improvement in accuracy over the base method for the same element number. The 
improvement is modest if the surface- vorticity method of generating a static 
solution (Figure 10b) is used, but is quite substantial if the interior suction 
method (Figure 10a) is used. 
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When the surface vorticlty solution for the Inlet of figure 11 was calcu- 
lated, the 244-element and the 149-element higher-order solutions gave graphically 
identical velocity distributions, which may be regarded as the ''correct" answer. 
This Is compared in Figure 13 with: the 103-element higher-order, the 103-element 

base method, and the 244-element base method. The 103-element higher-order 
solution is clearly superior to the 103-element base method and is at least 
comparable to and probably superior to the 244-element base method, which 
requires at least six times as much computing time. 

The suction method of simulating a static case (Figure 10a) was calculated 
with 149 elements for both the higher-order solution and the base method. For 
both solutions the suction element was placed at x = 44. The resulting veloc- 
ity distributions are compared with the "correct" (i.e., 244-element higher 
order) solution of Figure 14. It is evident that the higher-order solution 
effects a considerable increase In accuracy on the Inside of the inlet. In 
fact, the very close agreement of the two quite different higher-order solutions 
of Figure 14 is further evidence that the 244-element (or 149-element) higher- 
order solution with surface vorticity is indeed the "correct” static solution,, 
However, on the outside of the inlet the difficulties associated with suction 
solutions (Appendix D) are evident. Rather than a montonically decreasing vel- 
ocity of constant direction (towards the inlet lip) the suction solutions have 
stagnation points aft of which the flow is in the wrong direction. Moreover, 
the higher-order solution is not much better than the base method. Aft of the 
ficticious stagnation point of the static solution the calculated variation of 
surface pressure with mass flow ratio has the opposite sign from the true 
variation. Thus, if the outside of the inlet is of interest, the surface 
vorticity method of simulating a static solution should be used in preference 
to the suction method. The previous statement is true for all inlets where 
the outside surface is essentially parallel to the inside surface at aft loca- 
tions (Figure 11). It does not hold for "bellmouth" or "flush" inlets, which 
are characterized by the fact that the outside surface does not bend around 
and eventually assume a constant radius (constant value of y) but instead 
proceeds radially outward to large distances at a nearly constant value of x. 

For these inlets, the suction method and the surface vorticity method of simu- 
lating the static solution are essentially equivalent. 
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CONCLUSIONS 


The higher-order method (parabolic-element linear-source) is both faster 
and more accurate than the base method (straight-line-element constant-source) 
for exterior flow about simple bodies and for interior flow in ducts. The 
precise amount of improvement depends on the body shape and on the velocity 
component considered, but an order of magnitude improvement in both speed and 
accuracy has been obtained in apparently typical cases. 

For inlet flows the effectiveness of the higher-order method Is less 
dramatic, but it still appears to offer some Improvement over the base method. 

Use of a surface vorticity distribution to generate a static solution for 
Inlets is much preferable to use of an Interior suction surface. 

Satisfactory results can be obtained in'" inlet cases with fewer elements 
than is currently customary. 
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APPENDIX A 

ELEMENT GEOMETRY FORMULAS 



As in the present program the following geometric quantities associated 
with the j-th element are computed from two adjacent input points (x., yJ 
and (xj + i > y j+1 ). 

x o = i (x j + x j+l ) 
y o ' 2" ^ y j + y j+l ^ 

A =V ( Vl _X j )Z + (y j+l “ y j )Z (A-l ) 

X j+1 _X j 

C0Sa = 

A 

sina = ^j±LZi 

A 
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Now the element curvature must be computed. Half the value of curvature 
is stored.. The half curvature is denoted c* However, the actual curvature 
2c Is printed out on the preliminary output. This is the only addition to 
that output. The curvature of the jr-th element is computed as the geometric 
mean of the curvatures of two circles: the "left” circle through the points 

(xj„l » yj-l)* Cxj, yj), and (xu*|, yi+*|) and the "right* 1 circle through the 
points txj, yj), (xj+i, yj+l)» and Txj+o* yj+2)* If the "left" and "right" 
circles have curvatures of opposite sign, the element curvature is set equal to 
zero. The above scheme will not work for the first and last elements of a body. 
For the first element the curvature is set equal to that of the "right" circle 
and the (nonexistent) "left" circle is ignored. For the last element the curv- 
ature is set equal to that of the 'left" circle and the (nonexistent) "right" 
circle is ignored. This scheme is proper for the overwhelming majority of bodies: 
(1) finite body on axis, e.g., a sphere, (2) semi-infinite body on axis, e.g., 
a hemisphere cylinder, (3) a semi-infinite body off the axis, e.g., an inlet, 

(4) a duct, and (5) a ring airfoil with sharp trailing edge. The only exception 
is the smooth "donut type" body for which it would be preferrable to form a 
"left" circle for the first element using the second-to~last input point and a 
"right" circle for the last element using the second input point. However, this 
case does not seem to be important enough to justify an option,, 

The basic calculational unit of the above is the computation of the curv- 
ature of a circle through three given points. Let the points through which 
the circle goes be (x-j , y-|), (x2> y?)» (X3, y 3 ). The order of the points is 
important, because it determines the sign of tne curvature. In the present 
application the curvature is negative on convex portions of the body 0 The 
curvature 2c is 



(A-2) 


where 


D = 4[( x 2 -x 1 )(y 2 -y 3 ) - (x 2 -x 3 )(y 2 _ yp] 
R = ^(Dx 2 - D n ) 2 + (Dy 2 - D k ) 2 


(A-3) 


D n - 2[(x* + y^)(y 3 -y 2 ) + (x 2 + y 2 )(y 1 -y 3 ) + (x 3 + y 3 )(y 2 — y-,)] 
D k = 2[(x^ + y^)(x 2 _ x 3 ) + (x 2 + y 2 )(x 3 — x 1 ) + (x| + y|) (x 3 - x 2 )] 


Once the curvature is known the offset distance n 0 (see sketch above) is 
determined by passing a symmetric parabola with that curvature at the vertex 
through the endpoints. This gives 


n = -c 


(A-4) 
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Then the control point coordinates are 


v = y + n cos a 
J J o 0 

These now replace x 0 and y 0 , which are discarded* The control point (x, y) 
also serves as the origin of element coordi nates . 

The "half” arc length s 1 of the curved element is obtained from 

s’ = | (1 +1 cV) (A-6) 


The arc length as on the Basic Case output is 

as = 2s' (A— 7 ) 


All integrations versus arc length used to construct the induced velocity 
matrices have integration limits of -s' to +s‘. 

There is also an option to input curvatures for each element. If this 
option is used, the curvature computation is bypassed, and a table of curva- 
tures for all elements of all bodies is input. These are divided by two to 
produce the half curvatures c, which are stored with the basic geometric 
data for each element. Quantities saved for each element are 

x. y, c. A, sin a, cos a, s* (A-8) 
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APPENDIX B 

SIMQUUAR SUBELEMENT EXPANSIONS 

The logical structure of this calculation is described in [2], [6], [7]. 
First the half arc length s 1 is tested. If 

s' < 0.08 y (B-l ) 

i 

then the effect of the entire element is given by the singular subelement series 
(discussed below) and the argument of this series S' is 

S’ = sVy (B-2) 

If on the other hand, 

S' > 0.08 y 

then the effect of the "middle" of the element, i.e., the portion 

-0.08y < s < 0.08 y 

is given by the singular subelement series with argument 

S' = 0.08 

The "ends" of the element, i.e., the two portions 

-s’ < s < -0.08 y and 0.08 y < s < s* (B-6) 

are treated as two off-diagonal elements and their effects are computed by the 
numerical integration scheme of the section Induced Velocity Matrices. 

The required singular subelement series are logically similar to those of 
the base method, which they replace. These series are listed on the following 
pages, equations (B-9), (B-10), (B-ll), (B-12), (B-l 3) , (B-14), and (B-15). 

The only new parameter in these equations is 

C = cy (B-7) 

These equations lack the additional 2^ normal velocity term. Thus, these 
terms must be added to the coefficient or cj. (The coefficients of 1 J 
and ajv2) are not affected.) The terms to be added are: 


(B-3) 

(B-4) 

(B-5) 
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Ax i symmetric 


Cross Flow 


Vortex 


x : 2ir sin a 

y : -2n cos a 

x : 2tt sin a 

y : ~2 tt cos a 

e : 0 

x : 2 rr cos ot 

y : 2tt sin a 


(B~8) 


In the formulas below the notation of the left sides is chosen to be consistent 
with [6], [7], 

Axi symmetri c Source 



Cj | 2 sin a (cos a - 2cjs' + (y 2 Sln a 
+ T(T sin a 

4 cosot S' + ^ cosjajg sin 2 0 + 7 + 6 In I 


? 11 c ' 

COS a I sin a + —r + In 5- 
6 8 


sin 2 a + - 2 cos a + 8 In $' - In 8 - 4C(1 + cos a) 


)'• 




36 


2C 

T~ 


i 

3 “ 


2 2 

COS a — 2 sin a — 2C COS a 


,3 


+ crj^y^ | sin a (cos a — 2C j S’ 


(B-9) 
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ck | (2 sin 2 a + 2 In J— + 4C cos a) S' 

( l . 4 ? c * 

+ ( 2 sin a + 3 sin cx — 3 — 3 In g — 

+ Jg - cos a (6 sin 2 . + iy + 3 In 8) - 4 + 2 In S’ + 8C (1 - 2C cos a) )$ 


+ aj ]) y | 4 sin aS' + (sin a|-6 sin 2 . 


11 +6 In g — 


+ y C si na jcosa - C ) 


+ a^y 2 

0 J 


j f (s 1 n 2 a 


+ y + 1 n y- + 2C cos 


aj $ 1 


(B— 1 0) 


Cross-Flow Source 


("If) = a j | ^ s in a^cos a — 2C^S 1 

+ (2 sin. cos aj^2 sin 2 a — 9 - 6 In y-j 

c r 2 ? 

+ j 2 s1na (-6 sin a + 4 sin a — 5 + 6 In 8 — 4 cos a) 


+ 2 — 6 In S’ — 16C sin a (cos a— C) 
+ aP | 4 cos a S' + (|y cos a 


) s ' 


? <;» i 

-2 sin a + 3 + 6 In y~ 


+ |£ [cos 2 .- sin 2 a- 2C sin aj)s' 
+ af)/ | | sin. (cos a — 2Cj S ' 3 | 


(B-ll ) 
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oj ^2 sin 2 ct + 4 + 2 in |^- + 4c)s- 

+ (t 2 l sin2 “ ( 6 sin2 “ ~ 43 - 24 In |^-) + 9 + 27 In |^- 
+ g [ COS a(-3 sin 2 a + 7 + 5 In g— ) — 1 + 4C cos a(cos a — 2C) js |3 j 

+ oj^y I 4 sin as- + [-6 Sin 2 a + 29 + 30 In |^ 

+ J C sin a COS a — C js ,3 j 
+ o( 2) y Z | f (sin 2 a + | + In 1^- + 2C cos a)s' 3 j 

(y fe") p = j' 4 ! 1 + 1n |— ) s ' + (k[ 10 sir,2a + (6 sinZ “ _9) ln 1“ 

i COS a + 6 — 6 COS a In js’ 3 j 

j'Vj-f > 1 » -(f - 1 " r) s ' 3 J*«S 2) f z j-3-(l* r)s' 


(B-12) 


+ c 

9 


(B-13) 


Axi symmetric Vorticit.y 

(vortex) = f-2 sin 2 a +2 In g 4C cos ajS' 


( 


sin 2 a (-6 sin 2 a — 72 cos a + 23 + 12 In |— ) — g- (1 + In 


+i 


36 


, 2 


cos a(54 sin a - 1 1 + 6 In 8 + 30 In S' ) - 12 
+ 24C (3 sin 2 a - cos 2 a - C cos a) js' 3 (B-14) 
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(vortex) = ^2 sin a cos ce - 2cj)s' 

— sin a cos ct|-2 sin 2 a + 9 + 6 In |-j 

+ sin a|66 sin ot + 31 — 12 cos ct + ^ — In S* 

- 18 In 8 + 48C cos a]W 3 (B-15) 


24 



APPENDIX C 

NUMERICAL DIFFERENTIATION FORMULAS 

First the options must be applied as follows: 

Constant Source Option — set = 0 

- 42 ) 

Linear Source Option — set = 0 

■ J 


Now for any element that is not the first or last of a body, compute the 


quantities 







W here the "half arc length" s' of an element is defined in Appendix A. These 
are used in (11 ). 

When computing the effect of the first (j=l) element of each body, (11) 
is replaced by 

= Aa*| + Bc? 2 + Ca^ 

(C-2) 

°r = G 2 a l + H 2°2 + r 2°3 
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Here the quantities G 2 , H 2 , l? are given by (C-l) with j=2, i<,e os they are 
the same as the second elements (but they multiply different quantities). The 
new quantities are 


s' + + l/2(s‘ + s') 

A = " (Sg + s-pis!, + V2(s 1 r "+"s 3 TJ 


^ s£ + 1 /2 ( s-j + s|) 

8=2 TS\ + + ITJ 


(03) 


C = 2 


’(si + s 


3 1 

nr? 


+ s, 


+ 1 /2(s-j + 


spT 


When computing the effect of the last element on a body* say j=L, (11) 
is replaced by 


a l! 1) = Xa L-2 + Yct L- 1 + Zo L 
(2) 

°L “ G L-l a L-2 + H L-1 C L-1 + I L-l a L 


(C-4) 


where the quantities H^i* I|_-l are given by (C-l) with j=L-l, i.e., 
the same as the preceding element (but used differently). The new quantities 
are 


X = ^ 


K -2 + ^1 


L-l 

TR 


+ s, 


L-l 


l/Z(s -_ 2 + i[Tl 


„ S L-1 + 1/2 K-2 + S L ) 
Y = ’ 2 (s-_ 2 + iTTTT 5TT + 


’L-l ' ' L-l 




Z 


S L-1 + S L + 1/2(s L-l + S L-2 ) 
( S L-1 +s L)L' s L-l" V2 ( S L-2" +V L )J 


(C-5) 
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APPENDIX; D 

BASIC SOLUTIONS FOR INLETS AND SHROUDED PROPELLERS 

A common application of the present method is the calculation of flow over 
the forward portions of inlets or propeller shrouds. In such cases a complica- 
tion arises from the necessity of obtaining the desired total flow through the 
Inlet — the so-called mass flow ratio. (The effect of all the interior machin- 
ery is lumped Into this quantity.) The straightforward way to handle this prob- 
lem is to simply define by input points a surface across the interior coinciding 
with the actual propeller plane or compressor face and to specify the desired 
mass flow as a nonzero normal velocity on this surface. This procedure fails 
for reasons related to the Inability of the present method to handle internal 
corners. The calculated flow in the neighborhood of the simulated compressor 
face is erratic and not physically meaningful* Thus resort must be made to an 
artifice. 

The forward portion of the inlet is artificially extended by means of long 
afterbodies with constant inner and outer diameters as shown in figure 9. A 
forward location Is selected In the region where the flow is of interest, and 
it is designated the control station for mass flow ratio. Often the choice is 
the propeller plane or compressor face* Three fundamental flows are calculated. 
The first two are illustrated In figure 9. The long afterbody is left open at 
the rear and flows calculated for uniform onset flows at 0° and 90° angle of 
attack. The 0° flow gives a certain value of mass flow at the control station. 
This value cannot be predicted a priori, but it is usually near unity. The 
90° flow is a cross flow and gives zero net flow through the control station. 
These two solutions can be combined to give an onset flow at any angle of 
attack, but the mass flow ratio is always that obtained from the 0° onset flow. 
To obtain other mass flow ratios a third basic flow is required. 

The third basic flow is that for the inlet in static operation. That is, 
the flow is zero at infinity but has a finite mass flow at the control station. 
This flow can be linearly combined with the first two to give the flow about 
the inlet at any angle of attack and any mass flow ratio at the control station. 
Thus all solutions are obtained from the basic three. 

It is the third solution for static operation that leads to numerical dif- 
ficulty. Two methods of obtaining it are Illustrated in figure 10. The most 
straightforward procedure is that illustrated in figure 10a. A surface is 
placed across the interior of the inlet far back in the constant diameter region, 
and a nonzero normal velocity is prescribed on this "suction” surface. The 
flow In the neighborhood of the surface is meaningless as described above, but 
it smooths out upstream and is well-behaved at the control station. The mass 
flow obtained at the control station may be different than that specified on 
the suction surface, due to "leakage" caused by numerical errors. However, 
since this solution is to be linearly combined with other solutions, the exact 
value of mass flow is unimportant. This value is essentially absorbed in the 
combination constant. 
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The solution of figure IQa gives good results inside the inlet and around 
the lip, but inaccuracies enter in the exterior region. Because a finite 
number of surface elements are used, the exterior flow H $ees through" the 
inlet surface directly to the suction surface and is pulled towards that 
surface. Thus as shown in figure 10a there Is a stagnation point on the 
exterior surface. Fluid to the left enters the inlet mouth, while fluid on 
the right flows aft along the inlet. For the exact solution, which could be 
approached by using a very large number of elements, the Inlet surfaces "screen 
out" the suction element from the exterior flow, and the suction effect is felt 
solely through the inlet mouth. Thus the actual flow is to the left over the 
entire exterior surface with the velocity falling rapidly to zero with increase 
ing distance from the inlet mouth. The calculated velocity is not large, but 
It is in the wrong direction in this region. 

For the exterior region, the scheme shown in figure 10b gives a more 
accurate solution for static operation. A constant unit strength vorticity 
distribution is taken to lie on the surface. Of course, the normal velocity 
on the body due to the vorticity distribution is not zero, and a surface source 
distribution is also required. The flow field of the constant vorticity dis- 
tribution provides another onset flow to the body, and the values of source 
density on the elements are determined in the usual way to give zero normal 
velocity at all control points in the presence of the vorticity distribution. 
The only difference between this flow solution and the solution for a uniform 
onset flow Is the right side of the linear equations. (This device also pro- 
vides lift for airfoils [4], although for axi symmetric inlets the vorticity 
is ring vorticity.) It can be shown that if the afterbodies are infinitely 
long, even with finite-element lengths, the solution obtained in this manner 
is the correct static solution. Specifically, there is constant flow far back 
inside and outside a leftward flow that falls to zero with increasing distance 
from the inlet lip (figure 10b). The error in this case arises from the fact 
that long finite afterbodies lead to small nonzero velocities on the exterior 
surface, but this velocity is in the right direction. 
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Figure 2. An elementary arc. 
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Figure 4. Errors in surface velocity calculated by various methods for a sphere 
in a uniform onset flow parallel to the y-axis. (a) Velocity in the 
xy-plane. 
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Figure 6. Errors in calculated surface velocity on an 8-to-l prolate spheroid in a uniform onset flow 
parallel to the y-axis. (a) Velocity in the xy-plane. 
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Figure 8. Flow in a contracting duct of area ratio 16 with total velocity fluxes calculated at various 
locations by the higher-order and base method. 
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Angle of attack solutions for inlets, (a) 0° angle of attack, 
(b) 90° angle of attack. 
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Figure 10. Two methods for simulating flow about an inlet in static operation, 
(a) Interior suction, (b) Surface vorticity. 
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c) 103 ELEMENTS 

Figure 12. Three input point distributions on an inlet. 
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Figure 13. Calculated surface velocity distributions on an Inlet tn static operation simulated by a 
surface vortlclty solution. 
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Figure 14, Comparison of calculated surface velocity distributions on an inlet in static operation 
simulated by the two methods. 



